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We study the multi-species replicator model with linear fitness and random fitness ma- 
trices of various classes. By means of numerical resolution of the replicator equations, we 
determine the survival probability of a species in terms of its average interaction with 
the rest of the system. The role of the interaction pattern of the ecosystem in defining 
survival and extinction probabilities is emphasized. 
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1. Introduction 

Consider an ecosystem formed by TV biological species, with populations n\, . . . , njv- 
Over time scales where the effect of mutations can be neglected, the population of 
species i is governed by the evolution equation 



fm, (1.1) 



where /; is the (Fisherian) fitness of that species [|3[|8j]. Equation (1.1) assumes that 
evolution is continuous in time and that populations are large enough as to neglect 
the effect of random drift (9) . The interaction between species in the ecosystem is de- 
scribed by the fitnesses which generally depend on the whole set of populations. 
They may also include the effect of migration from and toward the system. It can 
be straightforwardly shown that the frequency of species i in the whole population, 
%i = rij/yij rij, satisfies the evolution equation 

X>£ I • (1.2) 

Assuming that the fitnesses can be given as functions of x%, . . . , Xn gives a closed 
set of equations for the frequencies, usually referred to as the replicator equations 
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Note that the replicator equations admit an alternative interpretation as the 
evolution equations for the populations -instead of the frequencies- of an TV-species 
ecosystem. Since for all times ^ i £j = 1 (see also Sect. ||), this would correspond to 
an ecosystem with fixed total population, presumably saturating the carrying ca- 
pacity of its environment. Within this interpretation, thus, the replicator equations 
describe an ensemble of biological species with global competition, which limits the 
growth of the whole population to the carrying capacity of the environment, and 
subject to interactions that affect the individual fitnesses. 

The replicator model has found many applications in biological evolutionary 
problems at various levels. It has been first considered in connection with the evo- 
lution of biomolecular concentrations in the early stages of life on the Earth [Q . In 
the frame of game theory, in close connection with biological evolution , it has 
been used to describe the dynamics of strategy frequencies in multi-player games 
[13 . As a basic model of an ecosystem, it was discussed by Schuster and Sigmund 
[12 , who coined the term replicator. The discrete-time version of Eqs. (~L2) [|ll| 
constitutes a well-known model for the evolution of genetic frequencies in asexual 
haploid populations || . Most of these applications assume that the fitnesses depend 
linearly on the frequencies, fi = J^ . that the replicator equations reduce 

to 

X! ai 3 X 3 ~ X! a 3kXjXk ) • (1-3) 

The matrix A = {a^ } is frequently called fitness matrix and, from the viewpoint 
of game theory, plays the role of a payoff matrix. It can be shown that the replica- 
tor equations with linear fitness are formally equivalent to a multi-species Lotka- 
Volterra model [^ . 

In spite of the large potential of the replicator equations to model evolutionary 
biological systems, the dynamical properties of their solutions are poorly under- 
stood. Besides the specific applications quoted above, generic mathematical prop- 
erties of the replicator model are known for very special situations only ]6|Jli|]. 
Typically, they involve a small number of biological species with quite unrealistic 
interaction patterns. Real ecosystems, on the other hand, usually involve tens of 
species with hundreds of food- web links [^J^] . Our aim in this paper is to approach 
this complex reality by exploring the replicator dynamics for large ecosystems with 
nontrivial interaction patterns. We deal with the problem from a statistical point 
of view and perform extensive numerical calculations, focusing on a probabilistic 
description of survival and extinction of species in terms of the parameters that 
define their interaction. 

In the next two sections we present -mainly for completeness and subsequent 
reference- various, rather general mathematical properties of the replicator equa- 
tions, including symmetry properties and stability conditions for special classes of 
random interaction patterns. Then, in the main part of the paper, we numerically 
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study more complex situations for the case of linear fitness, including multi-diagonal 
and disordered fitness matrices. Our results, which we interpret in connection with 
the average interaction of each species with the rest of the ecosystem, are summa- 
rized and discussed in the last section. 



2. Invariance properties and symmetries 

For future reference, let us begin discussing the invariance properties of the repli- 



cator model (1.2). In the first place, since the variables xi are defined as population 
fractions, Xi — rii/^-rij, the conditions Xi > and J2i x i = 1 must hold at 



all times. Equations (1.2) are compatible with these conditions, because the set 
Sn = {x — (xi, ...,Xfj)/xi > Vi and X4 = 1} is invariant under their action. 
In fact, the transverse component of x vanishes on the boundary of Sn, and 



i \ i / j 



(2.1) 



identically vanishes in Sn- Consequently, an initial condition in Sn generates a 
trajectory which remains within this set. From now on, we restrict the analysis of 



Eqs. (1.2) to Sn, where the dynamics is effectively N — 1-dimensional. In fact, Sn is 



an N — 1-dimcnsional polyhedron embedded in the ./V-dimcnsional Euclidean space. 



Equations (1.2) have moreover two important symmetries. First, multiplication 
of all the fitnesses by the same constant fi — > c\fi, is equivalent to a linear change 
of the temporal scale, t — > c\t. Second, addition of the same constant to all the 
fitnesses, fi — > fi + C2 leaves the equations invariant. This second symmetry implies 
that any meaningful result regarding the replicator dynamics will be given in terms 
of differences between fitnesses, rather than in terms of their individual values. 



Note, in fact, that Eqs. (1.2) can be rewritten as 

(2-2) 



For linear fitness, Eqs. (1.3), the above symmetries respectively amount to multi- 
plying and adding a constant to each element of the fitness matrix A, ctij — » Cidy 
and ctij — > ciij + C2 . 

3. Linear fitness: Ecological attitude and behavior 

In the following, we focus the analysis on the case of linear fitness, /, = a ij x j- 



For the sake of simplicity, the restriction of Eqs. (1.3) to Sn will be referred to as 
the replicator equations. 

The two summations in the right hand side of Eqs. (|l.3|) couple the dynamics 
of the species in different ways. The first summation, aijXj, coincides with the 
fitness of species i, and is a direct measure of the effect of the whole ecosystem 
on the reproductive success of that species. The coefficient a^, in fact, weights the 
contribution of species j to the fitness of i. The second summation, ^2 kl auXkXi is 
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the average fitness in the ecosystem, and acts as a global coupling. This "interac- 
tion" -a negative contribution for positive au~ is responsible for the preservation 
of the normalization of X^ 2 -*- The invariance of the system under summation of 
a constant to the fitness, quoted in Sect. || implies that all dynamical features in 
the replicator model depend on the fitness matrix A through differences between 
its elements. As a matter of fact, we find that a set of relevant quantities in the 
characterization of such features is given by the differences gij = a,ij — ajj. They 
compare the contribution of species j to the evolution of i with its self-contribution. 
We call cfcj the attitude of species j toward species i. A positive attitude, gij > 0, 
indicates that species j contributes to the growth of i more than to its own growth. 
It is therefore an altruistic attitude. On the other hand, g^ < stands for an 
egoistic attitude. By definition, the "self-attitudes" gjj are all zero. As we show in 



the following, the stability properties of many of the equilibria of Eqs. (1.3) can be 
given in terms of the attitudes gij . 



For large N, Eqs. (1.3) have generally a large number of fixed points. They can 
be classified by their number of non-vanishing coordinates, i.e. by the number of 
surviving species at each fixed point. There are N fixed points with only one surviv- 
ing species. The coordinates of these points are all zero except that corresponding 
to the surviving species, which equals unity. We emphasize that such equilibrium 
points, situated at the vertices of Sn, always exist, independently of the values of 
the coefficients a.y . In contrast, the number of fixed points with n surviving species 
(1 < n < N) depends on ay. We do know, however, that C(N, n) = N\/n\(N — n)\ 
is an upper bound for this number. These equilibrium points lie on the edges of 
the polyhedron Sn- Finally, only one equilibrium with N surviving species can 
exist [^). We call it the coexistence fixed point. If it exists, its coordinates satisfy 
J2j aijXj = J2m a kiXkXi for all i. 

For a general fitness matrix A and arbitrary N, the linear stability analysis 
of the fixed points of the replicator equations is quite a difficult task. Still, it is 
possible to find the general stability conditions for vertices. The j-th vertex, with 
Xj — 1 and Xi — for all i ^ j, is linearly stable if and only if g^ < for all i ^ j, 
i.e. if species j is egoistic toward all other species. Note that this condition can be 
simultaneously verified for several -or even all- species, in which case the system 
is multistable. 

An important case where the stability analysis can be fully achieved for all the 
fixed points is that of a diagonal matrix (a.^ = for i ^ j), which will be used as 
a reference case in our statistical analysis of survival and extinction. In this case, 
as above, the j-th vertex is linearly stable if and only if g^ < or, equivalently, 
ajj > 0. The equilibria on the edges of Sn, on the other hand, are always unstable. 
Meanwhile, the coexistence point is stable if and only if g^ — —o-jj > for all i, i.e. 
if all species are altruistic toward all other species. Consequently, if the coexistence 
point is stable, it is the only stable equilibrium, since the vertices are in turn 
unstable. The global stability of the coexistence point follows also from a general 
theorem for the replicator system with symmetric fitness matrices B. 
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Fig. 1. Dynamical behaviors for a two-species system: (a) coexistence, (b) bistability, and (c) 
domination by one species. In the three cases, empty dots indicate unstable fixed points, while 
full dots indicate stable equilibria. 

Linear stability analysis can also be carried out explicitly for small values of 
N. For two species, the dynamics is effectively one-dimensional. The phase space 
S2 can be mapped onto the interval [0,1], where one of the variables (say x\) 
takes on its values, while x-i = 1 — X\ at all times. In Fig. |l| we sketch the three 
possible dynamical scenarios. Coexistence (Fig. [|a) occurs when both species have 
positive attitude toward each other, while bistability (Fig. (TJa) occurs when both 
species have negative attitudes. Domination by one of the species takes place when 
one species is altruistic and the other species is egoistic. In such case the egoistic 
species dominates, while the other becomes extinct. For three species the dynamics 
is effectively bidimensional, and there is a wider spectrum of dynamical behaviors, 
including the possibility of existence of heteroclinic cycles . In general, we find 
that stability properties depend not only on the signs of the attitudes , but also 
on other combinations of the elements . 

Real ecosystems, however, are known to involve at least several species with 
complex patterns of interactions d|§. The restricting conditions that make the 
replicator model analytically solvable fall well apart from any situation with bio- 
logical significance. We therefore turn our attention to many-species systems with 
nontrivial interactions. In order to explore an ample class of fitness matrices, we 
choose the elements a,j at random, from prescribed distributions. This choice calls 
for a statistical analysis of the results, which in turn require resorting to the numer- 
ical resolution of the replicator equations. In the next section we discuss methods 
and present results of such analysis. 

4. Random fitness matrices for many-species ecosystems 

As advanced above, we analyze in this section the dynamics of the replicator model 
with a random fitness matrix and for a large number of species. In a generic sit- 
uation, the coexistence equilibrium will be unstable and a substantial number of 
species will become extinct. We are particularly interested in studying the depen- 
dence of the number of surviving species on the distribution of the elements a^-. 
Our ultimate goal would be to identify the ecological attitude that a species should 
have in order to optimize its probability of surviving within a given population, i. e. 
for a given distribution of . From then on, in order to make the description of 
interactions more compact, we characterize a species j by its mean attitude 

i i 
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Depending on the distribution of ay , the mean attitudes in a given population will 
take values in different ranges. For each realization of the evolution of the system, 
we respectively call n{G) and n s (G) the (normalized) distribution of the number 
of species and of the number of surviving species on G. Namely, for small 5G, 
the products Nn(G)8G and Nn s (G)5G represent, respectively, the total number 
of species and the number of surviving species in the interval (G, G + SG). The 
average fraction of surviving species as a function of the mean attitude is given by 

where (•) indicates average over realizations, with different ay and different initial 
conditions. The product p(G)5G measures the survival probability of a species with 
mean attitude between G and G + 5G within the given population. 



Taking into account the symmetries of the replicator equations (1.3), discussed 
in Sect. ^|, we have selected a set of distributions for ay that results to be repre- 
sentative and interesting. In most cases, we fix a significant subset of the matrix 
elements equal to zero, while the remaining ay are uniformly distributed on a finite 
interval. The subset of null elements can be fixed deterministically or at random 
over the matrix, and defines a reference level for all the other ay. Three intervals 
for the uniform distribution of the nonzero elements will be considered, namely, 
(— |, (~ \i |)i an d (■§, §)■ In this way, we take into account translations of the 
distribution of the nonzero elements with respect to the reference value. Note that, 
on the other hand, the three intervals have the same width. We briefly discuss the 



effect of changing this width at the end of Sect. 4.3. 

A special reference case which will also be considered in the following is a fit- 
ness matrix where no element is a priori fixed to zero, and ay is chosen from the 
uniform distribution for all i and j. This special case is statistically invariant under 
translations of the interval and changes of its width -up to a modification in the 
time scale. 



4.1. Numerical calculations 

In the numerical resolution of the replicator equations we face a main difficulty, 
namely, the dynamical instability of the invariant set Sjv- From Eqs. ( |2.1| ) it follows 
that, depending on the sign of the global coupling term J2j x jfj = k a jkXjXk, a 
small perturbation on the condition n = 1 may be amplified by the dynamics, 
driving the system away from Sn- This is expected to happen in the numerical 
calculations, due to rounding-off and discretization effects. Using a second-order 
Runge-Kutta (RK) integration algorithm |fTofl , we have observed departures from 
Sn for very long times only, typically when the evolution is very slow and the system 
is close to an unstable fixed point. However, such events prompt us to modify the 
integration algorithm, in order to ensure permanence on S n at all times. To solve 
the problem, we have included a normalization operation at each time step. After 
computing the updated variables Xi(t + dt) using the RK algorithm, we compute 
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the sum a = Xi(t+dt) and renormalize the variables as Xi(t + dt) — > Xi(t+dt)/cr. 
We have exhaustively tested this algorithm and observed that, in the regions where 
the original RK algorithm does not produce divergences, the solutions by both 
methods are indistinguishable. Furthermore, for the original RK algorithm we have 
observed that the perturbations that drive the solution away from Sn are initially 
of the order of the working precision (~ 10 -16 ), and the corrections introduced 
by the normalization procedure are of the same order, or ~ 10~ 15 at most. Such 
modifications are therefore much smaller than the discretization error of the RK 
algorithm, and the effect of normalization on the numerical precision results to be 
negligible. 

We have chosen initial conditions for x% with uniform distribution over Sjy. They 
have been generated by drawing N — 1 uniformly distributed random numbers j/j 
in (0, 1), with t/i < y 2 < ■ ■ ■ < 2/jv-i, and defining af»(0) = Vi+i - y x for 1 < i < N 
and x N (0) = 1 - y N -i- 

After integration over sufficiently long times, the values of X4 corresponding 
to the species which are in the way of becoming extinct have typically reached ex- 
tremely low levels. Such levels are usually meaningless from the biological viewpoint. 
In fact, the discreteness of biological populations fixes a lower bound of Q2j 71 j) _1 
for each Xi Q . Since the replicator model does not keep trace of the total population 
in the ecosystem, we choose as a threshold for extinctions x m \ n = 10~ 3 /iV. If, at 
the end of the integration, a variable Xi has reached a value below x m i n , species i is 
considered to have undergone extinction. To study the evolution of the number of 
surviving species, this definition is extended to all times. Surviving species at time 
t arc hence defined as those for which Xi(t) > x m i n . Our value of a; m i n is essentially 
arbitrary, but we find that the results are robust with respect to variations in this 
extinction threshold. 

Our calculations correspond mostly to the case N = 32, though we have verified 
that our main qualitative conclusions are valid for considerably larger ecosystems. 
As indicated above, our statistical analysis requires to average results over realiza- 
tions with different initial conditions and different matrix elements. Averages were 
taken over sets of 900 to 6000 realizations. In the following, we present the main 
results of our analysis, which focuses on the evolution of the number of surviving 
species and on the determination of the survival probability as a function of the 
mean attitude of each species. 

4.2. Diagonal and full-matrix systems 

We begin by analyzing, as a reference situation, the case of a diagonal fitness matrix 
(fly = for i ^ j). Recall that the stability properties in this case can be worked 
out exhaustively (Sect. ||). Note also that, for the diagonal system, = —ajj is 
independent of i, such that the attitude of each species j toward the remaining 
species is uniform. 

When the diagonal elements are drawn at random from then interval (— |, — i), 
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Fig. 2. (a) Number of surviving species as a function of time for diagonal fitness matrices and 
for the full-matrix case, with nonzero elements uniformly distributed on the indicated intervals, 
(b) Long time limit of the number of surviving species in the full-matrix case as a function of N. 



all these elements are negative. In this case all the species are altruistic, the coexis- 
tence point is globally stable, and all the species survive. The survival probability 
p(G) is thus a flat distribution. In the case that the diagonal elements are drawn 
from (— |, i), typically half of the species -those with a,jj < 0- are altruistic toward 
any other species, while the other half are egoistic. The only stable equilibria are the 
approximately N/2 vertices corresponding to the survival of a single species with 
cijj > 0. The surviving species is selected, in each realization, by the initial condi- 
tion. Finally, when the diagonal elements are chosen from (i, |), all the species are 
egoistic and the N vertices are stable. These are again the only stable equilibria and 
only one species survives. Figure ^a shows the evolution of the average number of 
surviving species for the above three cases. The corresponding survival probability 
distributions p(G) will be discussed later on, for comparison with other random 
matrices. 

A second case of reference is given, as advanced above, by the situation where 
no elements of the fitness matrix are a priori fixed to zero. We refer to this situation 
as the full-matrix case. The evolution of the average number of surviving species 
for the full-matrix case is also shown in Fig. 12^,. It is interesting to point out the 
relatively low number of species that ultimately survive in this case, which for 32 
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species is around 2.28 on the average. In Fig. ||b we show the average number of 
surviving species for the full-matrix system as a function of N. It does not change 
significantly with N, since for N — 128 it is around 2.35. As for the diagonal matrix, 
the distribution p(G) for the full-matrix case will be presented later. 

4.3. Systems with multi- diagonal matrices 

Taking the diagonal-matrix and the full-matrix systems as limiting reference cases, 
we now study intermediate, more complex situations. We consider Eqs. ( |l.3| ) for 
matrices with different numbers of nonzero diagonals. We define a fc-diagonal ran- 
dom matrix, with k an odd natural number, as a matrix whose elements a^- vanish 
for \i — j\ > (k — l)/2 and which are otherwise taken at random from a given dis- 
tribution. Clearly, k = 1 corresponds to the diagonal matrix, while the full-matrix 
case corresponds to k = 2N — 1. 

Figure ^a shows the evolution of the mean number of surviving species for multi- 
diagonal matrix systems with several values of k, and with the nonzero elements 
distributed uniformly on the interval (— |, — h). It can be seen that the average 
number of extinctions rapidly increases with the number of nonzero diagonals. In 
Fig. ^b we show the average normalized distribution of mean attitudes (n(G)) for 
three of the five cases analyzed in Fig. ||a. The range of mean attitudes where (n(G)) 
is sensibly different from zero moves toward the left as k is increased, reflecting the 
fact that the number of egoistic species grows. Combining this observation with the 
results of Fig. ||a we conclude that the more egoistic the population is, the less the 
number of surviving species. 

In Fig. |^c we show the numerical result for p(G) in the same realizations of 
Fig. ||a. Since in the diagonal-matrix case (k = 1) all the species survive, p(G) cor- 
responds to a uniform distribution. The results for k = 3 and k = 7 are remarkably 
different. The survival probabilities are in general lower than in the diagonal-matrix 
case and decrease almost linearly with the mean attitude over a wide range of values 
of G. Near the ends of the relevant intervals the linear behavior breaks down, and 
p(G) shows abrupt maxima. The maximum on the left side of the intervals indicate 
that, as expected, being one of the most egoistic species is a good strategy, as far as 
survival is concerned. Less obviously, the rightmost maximum indicates that being 
one of the most altruistic species is also convenient for survival. For k = 15, p{G) 
does not attains a maximum on the left extreme of the relevant interval and an 
intermediate optimum value of the mean attitude, close to G = appears, instead. 

In Figs. U and [5] we present the results for multi-diagonal matrices with their 
non- vanishing elements distributed uniformly on (— g, 5) and (4, §), respectively. In 
these situations, the population is more egoistic than in the case of Fig. ^, and the 
number of surviving species is very small. In contrast with that case, the number 
of survivals in the diagonal-matrix system is lower than for the full-matrix system. 
Now, only one species survives for a diagonal matrix, and the full-matrix case 
appears to be the situation with the highest rate of survivals. In the case of Fig. |J, 
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Fig. 3. (a) Number of surviving species as a function of time, (b) distributions of mean attitudes 
P(G), and (c) survival probability p(G), for multi-diagonal fitness matrices with nonzero elements 
taken from a homogeneous distribution on (— | , — ~). 



the transition between the limiting situations is smooth: the profiles of p(G) for 
different values of k are similar. All of them decrease monotonically with the mean 
attitude for intermediate values of G, and present abrupt maxima at the ends of 
the relevant intervals. The situation in Fig. ^|is more complex, and somehow more 
similar to the case of Fig. 0, since a maximum of p for an intermediate value of G 
is observed for relatively high values of k. 
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As a function of G, the fraction of survivals p characterizes the probability of 
surviving for a given species in terms of its ecological attitude toward other species 
in the system. It is also possible to characterize the survival probability in terms of 
the average attitude of the ecosystem toward a given species j, which we define as 

i i 

Thus, we have also studied pn(R), the average fraction of surviving species for a 
given value of R. Figure || presents the corresponding results for fc-diagonal systems 
with nonzero elements uniformly distributed on (— |, — |) and different values of k. 
We observe a monotonous growth of pr, clearly indicating that the probability of 
surviving for a given species is higher when the remaining species are more altruistic 
toward it. 

Up to now, we have kept the width of the distribution of nonzero elements 
equal to unity. Here, we briefly analyze the effect of changing this parameter in 
a particular situation. We consider a 7-diagonal matrix with its nonzero elements 
distributed uniformly in an interval of width A around — 1 . Numerical calculations 
show that the average number of extinctions as a function of time is practically 
independent of A. Figure shows p(G) for different values of A. As expected, 
the relevant interval of attitudes grows as the width of the distribution increases. 
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G 



Fig. 5. (a) Number of surviving species as a function of time and (b) survival probability p(G), 
for random multi-diagonal matrices with the nonzero elements uniformly distributed on ( ^ , § )• 
Note that, for 31 diagonals, more than 1/4 of the matrix elements are nonzero. Even though, the 
number of surviving species is nearly the same as in the full matrix case, which corresponds to 
k = 2N - 1 = 63. 

However, no qualitative changes are observed in the profile of p. 
4.4. Disordered random matrices 

Multi-diagonal fitness matrices represent interaction patterns for which the ensem- 
ble of species can be ordered in a well-defined sequence, where each species interacts 
with a few neighbors -besides the global interaction that couples the whole ecosys- 
tem. We now turn the attention to systems where the interaction pattern is more 
disordered than in the multi-diagonal case. We first consider fitness matrices where 
the diagonal elements are uniformly distributed in a certain interval, while the non- 
diagonal elements are drawn from the same distribution with probability p and 
fixed to zero with probability 1 — p. The diagonal-matrix case corresponds to p = 
and the full-matrix system is attained for p = 1. 

In Fig. H we show results for the average number of surviving species as a 
function of time, and for the distribution p(G), in a system where the nonzero 
matrix elements are taken from a uniform distribution on (— |, — A), and for different 
values of p. We find that the distributions p{G) for < p < 1 regularly interpolate 
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Fig. 6. Fraction of surviving species pn as function of R for a 7-diagonal matrix with nonzero 
elements uniformly distributed on(— i, ^). 
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Fig. 7. Survival probability p(G) corresponding to 7-diagonal matrices with nonzero elements 
uniformly distributed on intervals of different widths A around —1. 



the two extreme cases. The same feature is found when the nonzero elements are 
taken from (— |, |) and («, §). 

In Fig. |^ we compare the results for a 3-diagonal system and for a disordered 
system with probability p = 2/N = 1/16, which correspond to the same average 
number of nonzero matrix elements. It can be seen that the number of surviving 
species is significantly higher in the disordered system. Note that both systems 
have the same distribution of attitudes (n(G)). This is because, by definition, this 
function does not depend on the ordering of the interactions, but only on the average 
number, and distributions, of the nonzero matrix elements. Hence, the difference on 
the number of extinctions from both systems is a pure consequence of the disorder 
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of the interactions. 




Fig. 8. (a) Number of surviving species as a function of time and (b) survival probability p(G), 
for random matrices with diagonal elements taken from homogeneous distributions on (— i), 
and non-diagonal matrix elements taken from the same distribution with probability p and fixed 
as zero with probability 1 — p. Note that p = corresponds to the diagonal matrix case while 
p = 1 corresponds to the full matrix case. 



Finally, we study the replicator for a particular set of matrices that constitutes 
a continuous interpolation between multi-diagonal and disordered matrices. We 
consider a fc-diagonal fitness matrix and, with probability q, we set to zero each 
element of the k — 1 nonzero sub-diagonals. For each of these vanishing elements, 
we add a new random element outside the multi-diagonal block, \i—j\ > (k — l)/2. 
In this way it is possible, for instance, to transit continuously between the two cases 
reported in Fig. |^, starting from a 3-diagonal matrix (q = 0) and increasing the 
probability q to approach the disordered case (q = 1). Note that the distribution 
(n(G)) does not depend on q, and the only difference between systems with different 
values of this parameter is the degree of disorder in their interactions. In Fig. [n] we 
present results for this construction, but starting from a 7-diagonal matrix system. 
It can be seen that, as the degree of disorder in the fitness matrix grows, the number 
of survivals increases and the distributions p(G) become flatter. 
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Fig. 9. (a) Number of surviving species as a function of time and (b) survival probability p(G), for 
a 3-diagonal matrix whith nonzero elements uniformly distributed on (— | , — ^) and for a random 
disordered matrix with diagonal elements taken from the same distribution, and non-diagonal 
elements taken with probability p from the same distribution and fixed as zero with probability 
1-p. 



5. Conclusions 

In this paper we have presented a statistical study ol survival and extinction in 
the replicator model with linear fitness, exploring in particular the effect of the 
pattern of interactions in the ecosystem. We have selected a set of interaction pat- 
terns and fitness matrices that is expected to provide a representative (though 
not complete) sample of the different possible situations found in complex model 
ecosystems. Through the definition of the mean attitude of a species, we give a 
compact characterization of its ecological behavior, avoiding the problem of deal- 
ing with the detailed information involved in the fitness matrix. Our analysis has 
been focused on the study of the evolution of the number of surviving species for 
different interaction patterns, and on the determination of the survival probability 
as a function of the mean attitude. In the following, we summarize and comment 
our main conclusions. 

Concerning the number of surviving species in a given population, we have found 
that, in all cases, it increases as the mean attitude of the whole ensemble of species 
shifts to larger values. In other words, survival probabilities are higher in altruistic 
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Fig. 10. (a) Number of surviving species as a function of time and (b) survival probability p(G), 
for a 7-diagonal system that is being disordered by increasing the probability q. The nonzero 
matrix elements are uniformly distributed on (— §, — 5). 



populations. This property is clearly observed for multi-diagonal and disordered 
matrices, and for different distributions of the nonzero matrix elements. 

The survival probability grows also with the disorder of the fitness matrix, as 
illustrated by Figs. [l(] and |t], where we have compared the results for populations 
with different degrees of disorder in their interactions, but equal distributions of the 
mean attitude. In contrast with the case of multi-diagonal fitness matrices, where 
the number of interaction links is the same for all species, for disordered matrices 
there is a chance that the ecosystem becomes divided into weakly interacting sub- 
systems. Within one of such subsystems, each species interacts with a relatively 
smaller population and its survival probability thus increases. 

Once the fitness matrix of a population is fixed, we find that the survival prob- 
ability is higher for the most extreme attitudes. In most situations, in fact, species 
with the lowest and largest mean attitudes -i. e., the most egoistic and the most 
altruistic, respectively- have better chances of surviving. Some remarkable excep- 
tions have been found, however, where the egoistic extreme is not a good option, 
and an intermediate mean attitude maximizes the survival probability (see Fig. ^|). 
In contrast, the altruistic extreme has always been found to have a relative high 
chance, if not the highest, of surviving. 
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Let us finally stress that the present analysis has explored only a class of possible 
interaction patterns in the multi-species replicator model, with specific choices for 
the structure of the fitness matrix and the distribution of its elements. The possibil- 
ity is open for the consideration of much more general situations, approaching the 
statistical description of large real ecosystems with arbitrary, even highly complex, 
food-web structures. 
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